##Government Performance and Support for Democracy in Spain
#Last updated: Nov 10 2025

rm(list = ls())

##Load necessary packages for analysis:
library(foreign)
library(psych)
library(MASS)
library(ggplot2)
library(multilevel)
library(prediction)
library(dplyr)
library(car)
library(stargazer)
library(readr)
library(olsrr)
library(here)
library(rstudioapi)
library(tidyr)
library(colorspace)
library(GPArotation)
library(xtable)
library(margins)
library(ivreg)
library(dotwhisker)
library(readxl)
library(ggeffects)
library(purrr)
library(marginaleffects)
library(kableExtra)


################################
###To load extra functions for analysis:

rescale.to.1 <- function(var, num) (1 / (num - 1) * (var - num) + 1) #var=name of variable, num = max number on current scale (so, for a 7-point scale, you would put 7 -- starting number must be 1)


################################
##To load the data:

stub <- function() {}
thisPath <- function() {
  cmdArgs <- commandArgs(trailingOnly = FALSE)
  if (length(grep("^-f$", cmdArgs)) > 0) {
    # R console option
    normalizePath(dirname(cmdArgs[grep("^-f", cmdArgs) + 1]))[1]
  } else if (length(grep("^--file=", cmdArgs)) > 0) {
    # Rscript/R console option
    scriptPath <- normalizePath(dirname(sub(
      "^--file=",
      "",
      cmdArgs[grep("^--file=", cmdArgs)]
    )))[1]
  } else if (Sys.getenv("RSTUDIO") == "1") {
    if (rstudioapi::isAvailable(version_needed = NULL, child_ok = FALSE)) {
      # RStudio interactive
      dirname(rstudioapi::getActiveDocumentContext()$path)
    } else if (is.null(knitr::current_input(dir = TRUE)) == FALSE) {
      # Knit
      knitr::current_input(dir = TRUE)
    } else {
      # R markdown on RStudio
      getwd()
    }
  } else if (is.null(attr(stub, "srcref")) == FALSE) {
    # 'source'd via R console
    dirname(normalizePath(attr(attr(stub, "srcref"), "srcfile")$filename))
  } else {
    stop("Cannot find file path")
  }
}

nr.df <- read.csv(here::here(thisPath(), "Data","NatRep_12Aug22.csv"))
ov.df <- read.csv(here::here(thisPath(),"Data", "Oversample_12Aug22.csv"))

###############################################
##To merge, explore, and subset the data:
###############################################

#To merge the two datasets:
sp.df <- rbind(nr.df, ov.df)

##Subset the data to this analysis:
h.df <- subset(sp.df, H_mechanisms_1 > 0) #keep on only these rows
# Making 'h.df' the name of the dataframe for analysis

#look at columns
col_map <- data.frame(col = seq_along(nr.df), name = names(nr.df))
#View(col_map)

selected_columns <- c(1:47, 444:450, 478) #keeping only these columns
# Select columns using indexing
h.df <- h.df[, selected_columns]

##Control variables:
h.df$region <- car::recode(
  h.df$region,
  '"1"= "Andalucía"; 
                           "2"= "Aragon"; 
                           "3"= "Principado de Asturias"; 
                           "4"= "Illes Balears"; 
                           "5"= "Canarias"; 
                           "6"= "Cantabria"; 
                           "7"= "Castilla y León"; 
                           "8"= "Castilla-La Mancha"; 
                           "9"= "Catalunya"; 
                           "10"= "Comunitat Valenciana"; 
                           "11"= "Extremadura"; 
                           "12"= "Galicia"; 
                           "13"= "Madrid"; 
                           "14"= "Murcia"; 
                           "15"= "Navarra"; 
                           "16"= "País Vasco"; 
                           "17"= "La Rioja"; 
                           "18"= "Ceuta"; 
                           "19"= "Melilla"
                           ',
  as.factor = T
)

h.df$gender <- car::recode(
  h.df$gender,
  '
                           "1"= "Male"; 
                           "2"= "Female"
                           ',
  as.factor = T
)
h.df$age_recode <- car::recode(
  h.df$age,
  '
                               "4" = "1";
                               "5" = "2";
                               "6" = "3";
                               "7" = "4";
                               "8" = "5";
                               "9" = "6"
                               '
)
h.df$income <- h.df$SEL


##################
###To generate the treatment variable:
h.df$treat <- "Control"
h.df$treat[h.df$H_corruption != ""] <- "Corr Treat"
h.df$treat[h.df$H_policy != ""] <- "Pol Treat"
h.df$treat[h.df$H_corrupt_system != ""] <- "CS Treat"
h.df$treat[h.df$H_corrupt_elite != ""] <- "CE Treat"
h.df$treat <- as.factor(h.df$treat)
h.df$treat <- relevel(h.df$treat, ref = "Control")
table(h.df$treat)

##################
##Generate outcome variables:

##Democracy
h.df$H_support_dem_3_recode <- car::recode(
  h.df$H_support_dem_3,
  '
                                           "5"="1";
                                           "4"="2";
                                           "3"="3";
                                           "2"="4";
                                           "1"="5"
                                           '
)
h.df$H_support_dem_4_recode <- car::recode(
  h.df$H_support_dem_4,
  '
                                           "5"="1";
                                           "4"="2";
                                           "3"="3";
                                           "2"="4";
                                           "1"="5"
                                           '
)
h.df$H_support_dem_7_recode <- car::recode(
  h.df$H_support_dem_7,
  '
                                           "5"="1";
                                           "4"="2";
                                           "3"="3";
                                           "2"="4";
                                           "1"="5"
                                           '
)

##create index variable
h.df$dem16 <- (h.df$H_support_dem_1 + h.df$H_support_dem_6) / 2
hist(h.df$dem16)


##Recoding variables 0-1:
h.df$dem1 <- rescale.to.1(h.df$H_support_dem_1, 5)
h.df$dem2 <- rescale.to.1(h.df$H_support_dem_2, 5)
h.df$dem3 <- rescale.to.1(h.df$H_support_dem_3_recode, 5)
h.df$dem4 <- rescale.to.1(h.df$H_support_dem_4_recode, 5)
h.df$dem5 <- rescale.to.1(h.df$H_support_dem_5, 5)
h.df$dem6 <- rescale.to.1(h.df$H_support_dem_6, 5)
h.df$dem7 <- rescale.to.1(h.df$H_support_dem_7_recode, 5)
h.df$dem16 <- rescale.to.1(h.df$dem16, 5)


##Check how long it took people to complete the survey
describe(h.df$Duration..in.seconds. / 60)


##To remove speeders
h_speed.df <- h.df #dataset with speeders
hclean.df <- subset(h.df, (Duration..in.seconds. / 60) >= 4) #To drop speeders
h.df <- hclean.df #To run the analysis without speeders

# Create binary variables with NA for non-specific treatments
h.df$CE_Treat <- ifelse(
  h.df$treat == "CE Treat",
  "1",
  ifelse(h.df$treat == "Control", "0", NA)
)

h.df$Corr_Treat <- ifelse(
  h.df$treat == "Corr Treat",
  "1",
  ifelse(h.df$treat == "Control", "0", NA)
)

h.df$CS_Treat <- ifelse(
  h.df$treat == "CS Treat",
  "1",
  ifelse(h.df$treat == "Control", "0", NA)
)

h.df$Pol_Treat <- ifelse(
  h.df$treat == "Pol Treat",
  "1",
  ifelse(h.df$treat == "Control", "0", NA)
)

# Convert "0" and "1" to factor levels
h.df$CE_Treat <- factor(h.df$CE_Treat, levels = c("0", "1"))
h.df$Corr_Treat <- factor(h.df$Corr_Treat, levels = c("0", "1"))
h.df$CS_Treat <- factor(h.df$CS_Treat, levels = c("0", "1"))
h.df$Pol_Treat <- factor(h.df$Pol_Treat, levels = c("0", "1"))

##To merge the compliance coding back into the dataframe:
compliance.df <- read.csv(here::here(thisPath(),"Data","ComplianceCheckTextsR3(Sheet 1 - ComplianceCheckTextsR).csv"))
compliance.df$ResponseId <- compliance.df$ResponseID
h2.df <- merge(h.df,compliance.df,by="ResponseId",all.y=T)

################################
##Basic Analysis Models:
################################
#Dem16
summary(
  dem16.lm <- lm(dem16 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 1
summary(
  dem1.lm <- lm(dem1 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 2
summary(
  dem2.lm <- lm(dem2 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 3
summary(
  dem3.lm <- lm(dem3 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 4
summary(
  dem4.lm <- lm(dem4 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 5
summary(
  dem5.lm <- lm(dem5 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 6
summary(
  dem6.lm <- lm(dem6 ~ treat + region + gender + age + income, data = h.df)
)

##democracy 7
summary(
  dem7.lm <- lm(dem7 ~ treat + region + gender + age + income, data = h.df)
) 

#######################################
#####TO GENERATE FIGURE 1##############
#######################################

h_control <- subset(h.df, treat == "Control")
stacked <- gather(
  h_control,
  key = measure,
  value = Rate,
  c(dem16, dem2, dem3, dem4, dem5, dem7)
)
stacked$Rate <- as.numeric(stacked$Rate)
table(stacked$Rate)

stacked$agree_name <- ""
stacked$agree_name[stacked$Rate == 0] <- "Strongly Disagree"
stacked$agree_name[stacked$Rate == 0.125] <- "Strongly Disagree"
stacked$agree_name[stacked$Rate == 0.25] <- "Somewhat Disagree"
stacked$agree_name[stacked$Rate == 0.375] <- "Somewhat Disagree"
stacked$agree_name[stacked$Rate == 0.5] <- "Neither Agree nor Disagree"
stacked$agree_name[stacked$Rate == 0.625] <- "Somewhat Agree"
stacked$agree_name[stacked$Rate == 0.75] <- "Somewhat Agree"
stacked$agree_name[stacked$Rate == 0.875] <- "Strongly Agree"
stacked$agree_name[stacked$Rate == 1] <- "Strongly Agree"

#Generating plot for figure 1
ggplot(data = stacked) +
  geom_bar(
    aes(
      measure,
      fill = factor(
        agree_name,
        levels = c(
          "Strongly Agree",
          "Somewhat Agree",
          "Neither Agree nor Disagree",
          "Somewhat Disagree",
          "Strongly Disagree"
        )
      )
    ),
    position = "fill"
  ) +
  labs(fill = "", x = "", y = "", title = "") +
  scale_x_discrete(
    labels = c(
      "Conceptual Democracy",
      "Freedom of Assembly",
      "Freedom of Press",
      "Judicial Independence",
      "Rule of Law",
      "Party Freedom"
    ),
    guide = guide_axis(angle = 30)
  ) +
  scale_y_continuous(labels = c("0%", "25%", "50%", "75%", "100%")) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18),
    plot.subtitle = element_text(hjust = 0.5, size = 16),
    axis.text.x = element_text(size = 11),
    axis.title.x = element_text(size = 14),
    axis.text.y = element_text(size = 14),
    axis.title.y = element_text(size = 16),
    legend.text = element_text(size = 14),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.7)
  ) +
  scale_fill_discrete_diverging(palette = "Cork")
ggsave(
  here::here(thisPath(), "Figures", "Figure1.pdf"),
  width = 9,
  height = 6,
  family = "Times"
)
ggsave(here::here(thisPath(), "Figures", "Figure1.jpeg"), width = 9, height = 6)


#######################################
#####TO GENERATE FIGURE 2##############
#######################################

fig1 <- dwplot(
  list(dem7.lm, dem5.lm, dem4.lm, dem3.lm, dem2.lm, dem16.lm),
  dot_args = list(size = 3),
  whisker_args = list(size = 1)
) +
  scale_y_discrete(
    limit = c(
      "treatCE Treat",
      "treatCS Treat",
      "treatPol Treat",
      "treatCorr Treat"
    ),
    labels = c( 
      "Elite Corruption",
      "Systemic Corruption",
      "Unemployment",
      "General Corruption"
    )
  ) +
  xlim(c(-0.15, 0.15)) +
  labs(
    x = "Marginal Change in Support for Democracy",
    title = "",
    y = "Treatment"
  ) +
  geom_vline(xintercept = 0, linetype = 2, colour = "#234230") +
  scale_color_manual(
    values = c("#265386", "#7B9EDB", "#93A7BF", "#A7BFA8", "#4C8C67", "#234230")
  ) +
  annotate(
    "text",
    x = -0.014,
    y = 4.18,
    label = "Party Freedom",
    col = "#234230",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -0.04,
    y = 4.11,
    label = "Rule of Law",
    col = "#4C8C67",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -0.06,
    y = 4.045,
    label = "Judicial Independence",
    col = "#A7BFA8",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -0.013,
    y = 3.975,
    label = "Freedom of Press",
    col = "#93A7BF",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -0.057,
    y = 3.911,
    label = "Freedom of Assembly",
    col = "#7B9EDB",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = -0.07,
    y = 3.843,
    label = "Conceptual Democracy",
    col = "#265386",
    size = 3,
    hjust = 1
  ) +
  theme_bw() +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    axis.text = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 12),
    axis.title.y = element_text(size = 12)
  )
fig1 
ggsave(
  here::here(thisPath(), "Figures", "Figure2.pdf"),
  width = 8,
  height = 8,
  family = "Times"
)
ggsave(here::here(thisPath(), "Figures", "Figure2.jpeg"), width = 8, height = 8)


###############################
####Appendix Tables and Figures############
###############################

##Code to generate Appendix Table 1
table(h.df$region)
sum(table(h.df$region))

table(h.df$gender)
sum(table(h.df$gender))

table(h.df$age)
sum(table(h.df$age))

table(h.df$income)
sum(table(h.df$income))


#########################
##Code to generate Appendix Figure 1 and Table 2
demindex.df <- data.frame(
  h.df$H_support_dem_1,
  h.df$H_support_dem_2,
  h.df$H_support_dem_3_recode,
  h.df$H_support_dem_4_recode,
  h.df$H_support_dem_5,
  h.df$H_support_dem_6,
  h.df$H_support_dem_7_recode
)
##Then ran the simple function with the new dataframe.
cronbach(demindex.df) #.52

#Figure 1
vss_demindex <- vss(
  demindex.df,
  fm = "minres",
  rotate = "oblimin",
  title = "VSS: Demindex"
)
vss_demindex
pdf(
  here::here(thisPath(), "Figures", "Appendix_Figure1.pdf"),
  width = 8,
  height = 5,
  family = "Times"
)
VSS.plot(vss_demindex, title = "Factor Analysis: Democracy Index")
dev.off()

#Table 2: EFA results
demindex.mr_b <- fa(
  demindex.df,
  2,
  fm = "minres",
  rotate = "oblimin",
  scores = "Bartlett"
)
demindex.mr_b
names = c(
  "Democracy Better",
  "Freedom of Assembly",
  "Freedom of Press",
  "Judicial Independence",
  "Rule of Law",
  "Institutional Processes",
  "Pol. Parties"
)
table = cbind(names, round(as.numeric(demindex.mr_b$loadings), digits = 4))
table = as.data.frame(table)
table$Item = table$names
table$Loadings = table$V2
table$names = NULL
table$V2 = NULL
row.names(table) <- NULL
table = print(table, row.names = F)
print(xtable(table), include.rownames = F)


###############################
#Appendix A3- Main analysis table without speeders
stargazer(
  dem16.lm,
  dem2.lm,
  dem3.lm,
  dem4.lm,
  dem5.lm,
  dem7.lm,
  covariate.labels = c(
    "Intercept",
    "CE Treatment",
    "Corrupt Treatment",
    "CS Treatment",
    "Pol Treatment",
    "Aragón",
    "Principado de Asturias",
    "Illes Balears",
    "Canarias",
    "Cantabria",
    "Castilla y León",
    "Castilla-La Mancha",
    "Catalunya",
    "Comunitat Valenciana",
    "Extremadura",
    "Galicia",
    "Madrid",
    "Murcia",
    "Navarra",
    "País Vasco",
    "La Rioja",
    "Ceuta",
    "Melilla",
    "Gender (male)",
    "age",
    "income"
  ),
  type = "html",
  style = "apsr",
  title = "Support for Democracy",
  intercept.bottom = FALSE,
  intercept.top = TRUE,
  digits = 3,
  star.cutoffs = c(.05, .01, .001),
  dep.var.labels = c(
    "Institutional Democracy",
    "Freedom of Assembly",
    "Freedom of Press",
    "Judicial Independence",
    "Rule of Law",
    "Freedom of Political Parties"
  ),
  notes = c("Base categories: Control, Andalucia, Socialist Party, Catholic"),
  out = (here::here(thisPath(), "Tables", "Appendix_Table3.doc"))
)


###############################
#Appendix A4- Main analysis table with speeders

##Analysis with the speeders
#Dem16
summary(
  dem16s.lm <- lm(
    dem16 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 1
summary(
  dem1s.lm <- lm(
    dem1 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 2
summary(
  dem2s.lm <- lm(
    dem2 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 3
summary(
  dem3s.lm <- lm(
    dem3 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 4
summary(
  dem4s.lm <- lm(
    dem4 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 5
summary(
  dem5s.lm <- lm(
    dem5 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 6
summary(
  dem6s.lm <- lm(
    dem6 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

##democracy 7
summary(
  dem7s.lm <- lm(
    dem7 ~ treat + region + gender + age + income,
    data = h_speed.df
  )
)

#Creating full table
stargazer(
  dem16s.lm,
  dem2s.lm,
  dem3s.lm,
  dem4s.lm,
  dem5s.lm,
  dem7s.lm,
  covariate.labels = c(
    "Intercept",
    "CE Treatment",
    "Corrupt Treatment",
    "CS Treatment",
    "Pol Treatment",
    "Aragón",
    "Principado de Asturias",
    "Illes Balears",
    "Canarias",
    "Cantabria",
    "Castilla y León",
    "Castilla-La Mancha",
    "Catalunya",
    "Comunitat Valenciana",
    "Extremadura",
    "Galicia",
    "Madrid",
    "Murcia",
    "Navarra",
    "País Vasco",
    "La Rioja",
    "Ceuta",
    "Melilla",
    "Gender (male)",
    "age",
    "income"
  ),
  type = "html",
  style = "apsr",
  title = "Support for Democracy",
  intercept.bottom = FALSE,
  intercept.top = TRUE,
  digits = 3,
  star.cutoffs = c(.05, .01, .001),
  dep.var.labels = c(
    "Institutional Democracy",
    "Freedom of Assembly",
    "Freedom of Press",
    "Judicial Independence",
    "Rule of Law",
    "Freedom of Political Parties"
  ),
  notes = c("Base categories: Control, Andalucia, Socialist Party, Catholic"),
  out = (here::here(thisPath(), "Tables", "Appendix_Table4.doc"))
)


#############################
###To estimate and plot CACE effects -- Appendix Figure A.2:
#############################

################
####NEGATIVE####
################
## We want the control condition to be zero:
h2.df$Negative[is.na(h2.df$Negative)] <- 0
h2.df$On_topic[is.na(h2.df$On_topic)] <- 0
h2.df$Political_performance[is.na(h2.df$Political_performance)] <- 0

h2.df$Compliant <- 0
h2.df$Compliant[h2.df$Negative == 1 & h2.df$On_topic == 1] <- 1

table(h2.df$Negative)
table(h2.df$On_topic)
table(h2.df$Political_performance)
table(h2.df$Compliant)


#######################
##To estimate CE treat treatment effects for full compliers:
summary(
  cace_ce_dem16.lm <- ivreg(
    dem16 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem2.lm <- ivreg(
    dem2 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem3.lm <- ivreg(
    dem3 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem4.lm <- ivreg(
    dem4 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem5.lm <- ivreg(
    dem5 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem6.lm <- ivreg(
    dem6 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_ce_dem7.lm <- ivreg(
    dem7 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CE_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

#Estimating marginal effects
cace_ce_dem16.df <- margins_summary(
  cace_ce_dem16.lm,
  data = cace_ce_dem16.lm$model
)[4, ]
cace_ce_dem2.df <- margins_summary(
  cace_ce_dem2.lm,
  data = cace_ce_dem2.lm$model
)[4, ]
cace_ce_dem3.df <- margins_summary(
  cace_ce_dem3.lm,
  data = cace_ce_dem3.lm$model
)[4, ]
cace_ce_dem4.df <- margins_summary(
  cace_ce_dem4.lm,
  data = cace_ce_dem4.lm$model
)[4, ]
cace_ce_dem5.df <- margins_summary(
  cace_ce_dem5.lm,
  data = cace_ce_dem5.lm$model
)[4, ]
cace_ce_dem7.df <- margins_summary(
  cace_ce_dem7.lm,
  data = cace_ce_dem7.lm$model
)[4, ]


#######################
##To estimate CS treat treatment effects for compliers:
summary(
  cace_cs_dem16.lm <- ivreg(
    dem16 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_cs_dem2.lm <- ivreg(
    dem2 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_cs_dem3.lm <- ivreg(
    dem3 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_cs_dem4.lm <- ivreg(
    dem4 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_cs_dem5.lm <- ivreg(
    dem5 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)


summary(
  cace_cs_dem7.lm <- ivreg(
    dem7 ~
      On_topic +
      region +
      gender +
      age +
      income |
      CS_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)
#Estimating marginal effects for CS treat
cace_cs_dem16.df <- margins_summary(
  cace_cs_dem16.lm,
  data = cace_cs_dem16.lm$model
)[4, ]
cace_cs_dem2.df <- margins_summary(
  cace_cs_dem2.lm,
  data = cace_cs_dem2.lm$model
)[4, ]
cace_cs_dem3.df <- margins_summary(
  cace_cs_dem3.lm,
  data = cace_cs_dem3.lm$model
)[4, ]
cace_cs_dem4.df <- margins_summary(
  cace_cs_dem4.lm,
  data = cace_cs_dem4.lm$model
)[4, ]
cace_cs_dem5.df <- margins_summary(
  cace_cs_dem5.lm,
  data = cace_cs_dem5.lm$model
)[4, ]
cace_cs_dem7.df <- margins_summary(
  cace_cs_dem7.lm,
  data = cace_cs_dem7.lm$model
)[4, ]


#######################
##To estimate Corr treat treatment effects for compliers:
summary(
  cace_corr_dem16.lm <- ivreg(
    dem16 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_corr_dem2.lm <- ivreg(
    dem2 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_corr_dem3.lm <- ivreg(
    dem3 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_corr_dem4.lm <- ivreg(
    dem4 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_corr_dem5.lm <- ivreg(
    dem5 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)


summary(
  cace_corr_dem7.lm <- ivreg(
    dem7 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Corr_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

#Estimating marginal effects for Corr treat
cace_corr_dem16.df <- margins_summary(
  cace_corr_dem16.lm,
  data = cace_corr_dem16.lm$model
)[4, ]
cace_corr_dem2.df <- margins_summary(
  cace_corr_dem2.lm,
  data = cace_corr_dem2.lm$model
)[4, ]
cace_corr_dem3.df <- margins_summary(
  cace_corr_dem3.lm,
  data = cace_corr_dem3.lm$model
)[4, ]
cace_corr_dem4.df <- margins_summary(
  cace_corr_dem4.lm,
  data = cace_corr_dem4.lm$model
)[4, ]
cace_corr_dem5.df <- margins_summary(
  cace_corr_dem5.lm,
  data = cace_corr_dem5.lm$model
)[4, ]
cace_corr_dem7.df <- margins_summary(
  cace_corr_dem7.lm,
  data = cace_corr_dem7.lm$model
)[4, ]


##To estimate Pol treat treatment effects for compliers:
summary(
  cace_pol_dem16.lm <- ivreg(
    dem16 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_pol_dem2.lm <- ivreg(
    dem2 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_pol_dem3.lm <- ivreg(
    dem3 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_pol_dem4.lm <- ivreg(
    dem4 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

summary(
  cace_pol_dem5.lm <- ivreg(
    dem5 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)


summary(
  cace_pol_dem7.lm <- ivreg(
    dem7 ~
      On_topic +
      region +
      gender +
      age +
      income |
      Pol_Treat + region + gender + age + income,
    data = h2.df
  ),
  diagnostics = T
)

#Estimating marginal effects for Pol treat
cace_pol_dem16.df <- margins_summary(
  cace_pol_dem16.lm,
  data = cace_pol_dem16.lm$model
)[4, ]
cace_pol_dem2.df <- margins_summary(
  cace_pol_dem2.lm,
  data = cace_pol_dem2.lm$model
)[4, ]
cace_pol_dem3.df <- margins_summary(
  cace_pol_dem3.lm,
  data = cace_pol_dem3.lm$model
)[4, ]
cace_pol_dem4.df <- margins_summary(
  cace_pol_dem4.lm,
  data = cace_pol_dem4.lm$model
)[4, ]
cace_pol_dem5.df <- margins_summary(
  cace_pol_dem5.lm,
  data = cace_pol_dem5.lm$model
)[4, ]
cace_pol_dem7.df <- margins_summary(
  cace_pol_dem7.lm,
  data = cace_pol_dem7.lm$model
)[4, ]

# Combine all CACE results into a single data frame
cace_combined.df <- data.frame(bind_rows(
  cace_ce_dem16.df %>% mutate(Treatment = "CE Treat", Outcome = "dem16"),
  cace_ce_dem2.df %>% mutate(Treatment = "CE Treat", Outcome = "dem2"),
  cace_ce_dem3.df %>% mutate(Treatment = "CE Treat", Outcome = "dem3"),
  cace_ce_dem4.df %>% mutate(Treatment = "CE Treat", Outcome = "dem4"),
  cace_ce_dem5.df %>% mutate(Treatment = "CE Treat", Outcome = "dem5"),
  cace_ce_dem7.df %>% mutate(Treatment = "CE Treat", Outcome = "dem7"),

  cace_cs_dem16.df %>% mutate(Treatment = "CS Treat", Outcome = "dem16"),
  cace_cs_dem2.df %>% mutate(Treatment = "CS Treat", Outcome = "dem2"),
  cace_cs_dem3.df %>% mutate(Treatment = "CS Treat", Outcome = "dem3"),
  cace_cs_dem4.df %>% mutate(Treatment = "CS Treat", Outcome = "dem4"),
  cace_cs_dem5.df %>% mutate(Treatment = "CS Treat", Outcome = "dem5"),
  cace_cs_dem7.df %>% mutate(Treatment = "CS Treat", Outcome = "dem7"),

  cace_corr_dem16.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem16"),
  cace_corr_dem2.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem2"),
  cace_corr_dem3.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem3"),
  cace_corr_dem4.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem4"),
  cace_corr_dem5.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem5"),
  cace_corr_dem7.df %>% mutate(Treatment = "Corr Treat", Outcome = "dem7"),

  cace_pol_dem16.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem16"),
  cace_pol_dem2.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem2"),
  cace_pol_dem3.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem3"),
  cace_pol_dem4.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem4"),
  cace_pol_dem5.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem5"),
  cace_pol_dem7.df %>% mutate(Treatment = "Pol Treat", Outcome = "dem7")
))

# Add significance stars based on p-values
cace_combined.df <- cace_combined.df %>%
  mutate(
    star = case_when(
      #p <= .1 & p > .05 ~ "+",
      p <= .05 & p > .01 ~ "*",
      p <= .01 & p > .001 ~ "**",
      p <= .001 ~ "***"
    )
  )

cace_combined.df$Treatment <- car::recode(
  cace_combined.df$Treatment,
  '
                           "CE Treat"= "Elite Corruption";
                           "CS Treat"= "Systemic Corruption";
                           "Corr Treat"= "General Corruption";
                           "Pol Treat"= "Unemployment"
                           ',
  as.factor = T
)

cace_combined.df$Treatment <- factor(
  cace_combined.df$Treatment,
  levels = c(
    "Elite Corruption",
    "Systemic Corruption",
    "Unemployment",
    "General Corruption"
  )
)

#######################################
# Create Appendix Figure A.2: CACE Plot
#######################################
cace_plot <- ggplot(
  data = cace_combined.df,
  aes(x = Treatment, y = AME, color = Outcome, shape = Outcome, group = Outcome)
) +
  geom_point(
    stat = "identity",
    position = position_dodge(width = .6, preserve = "total")
  ) +
  geom_linerange(
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = .6, preserve = "total"),
    linewidth = 1
  ) +
  #geom_text(data = cace_combined.df, aes(label = star, group = Outcome), position = position_dodge2(width = .5, preserve = "total"), show.legend = FALSE) +
  theme(text = element_text(family = "serif")) +
  coord_flip() +
  theme(axis.title = element_text(size = 8)) +
  theme(axis.text = element_text(size = 8)) +
  #theme(legend.text = element_text(size = 7)) +
  #theme(legend.title = element_text(size = 7)) +
  geom_hline(yintercept = 0, lty = 2) +
  theme_bw() +
  theme(legend.position = "none", panel.grid.minor.x = element_blank()) +
  scale_color_manual(
    values = c("#265386", "#7B9EDB", "#93A7BF", "#A7BFA8", "#4C8C67", "#234230")
  ) +
  labs(y = "CACE: Marginal Change in Support for Democracy", x = "Treatment") +
  annotate(
    "text",
    x = 4.255,
    y = -.01,
    label = "Party Freedom",
    col = "#234230",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = 4.163,
    y = -.04,
    label = "Rule of Law",
    col = "#4C8C67",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = 4.058,
    y = -.0635,
    label = "Judicial Independence",
    col = "#A7BFA8",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = 3.96,
    y = -.014,
    label = "Freedom of Press",
    col = "#93A7BF",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = 3.858,
    y = -.0635,
    label = "Freedom of Assembly",
    col = "#7B9EDB",
    size = 3,
    hjust = 1
  ) +
  annotate(
    "text",
    x = 3.758,
    y = -.0765,
    label = "Conceptual Democracy",
    col = "#265386",
    size = 3,
    hjust = 1
  ) +
  scale_y_continuous(breaks = c(-.1, -.05, 0, .05, .1))
# Print and save the plot
cace_plot
ggsave(
  here::here(thisPath(), "Figures", "Appendix_Figure2_CACE.pdf"),
  width = 8,
  height = 8,
  family = "Times"
)
ggsave(
  here::here(thisPath(), "Figures", "Appendix_Figure2_CACE.jpeg"),
  width = 8,
  height = 8
)

#############################################
##### HETEROGENEOUS EFFECTS BY IDEOLOGY (Appendix Table 7)
#############################################

table(h.df$ideology_1, useNA = "ifany")

## Create left, right, middle variable
h.df$ideo_left_right_middle <- case_when(
  h.df$ideology_1 %in% c(0, 1) ~ "Left",
  h.df$ideology_1 %in% c(9, 10) ~ "Right",
  h.df$ideology_1 %in% 2:8 ~ "Center",
  TRUE ~ NA_character_
)

table(h.df$ideo_left_right_middle, useNA = "ifany")

outcomes <- c("dem16", "dem2", "dem3", "dem4", "dem5", "dem7")

###############################
## Interaction models for left, right, middle (Appendix Table 7)
###############################

# Fit interaction models for ideology
summary(
  dem16_left.lm <- lm(
    dem16 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem2_left.lm <- lm(
    dem2 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem3_left.lm <- lm(
    dem3 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem4_left.lm <- lm(
    dem4 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem5_left.lm <- lm(
    dem5 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem7_left.lm <- lm(
    dem7 ~ treat * ideo_left_right_middle + region + gender + age + income,
    data = h.df
  )
)

# Calculate average marginal effects of treatment by ideology
mfx_dem16 <- avg_slopes(
  dem16_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)
mfx_dem2 <- avg_slopes(
  dem2_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)
mfx_dem3 <- avg_slopes(
  dem3_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)
mfx_dem4 <- avg_slopes(
  dem4_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)
mfx_dem5 <- avg_slopes(
  dem5_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)
mfx_dem7 <- avg_slopes(
  dem7_left.lm,
  variables = "treat",
  by = "ideo_left_right_middle"
)

print(mfx_dem16)

# Combine into a single table
all_mfx <- rbind(
  mfx_dem16 %>% mutate(outcome = "Conceptual Democracy"),
  mfx_dem2 %>% mutate(outcome = "Freedom of Assembly"),
  mfx_dem3 %>% mutate(outcome = "Freedom of Press"),
  mfx_dem4 %>% mutate(outcome = "Judicial Independence"),
  mfx_dem5 %>% mutate(outcome = "Rule of Law"),
  mfx_dem7 %>% mutate(outcome = "Party Freedom")
)


# Clean up prime names for table
all_mfx <- all_mfx %>%
  mutate(
    contrast = case_when(
      grepl("CE Treat", contrast) ~ "Elite Corruption",
      grepl("Corr Treat", contrast) ~ "General Corruption",
      grepl("CS Treat", contrast) ~ "Systemic Corruption",
      grepl("Pol Treat", contrast) ~ "Unemployment",
      TRUE ~ contrast
    )
  )

# Format table for HTML output
formatted_table <- all_mfx %>%
  mutate(
    stars = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ ""
    ),
    estimate_display = sprintf("%.3f%s", estimate, stars),
    se_display = sprintf("(%.3f)", std.error),
    p_display = sprintf("%.3f", p.value)
  ) %>%
  select(
    outcome,
    ideo_left_right_middle,
    contrast,
    estimate_display,
    se_display,
    p_display
  )

#######################################
#####TO GENERATE Appendix Table 7
#######################################
html_table_ideo <- formatted_table %>%
  kable(
    format = "html",
    col.names = c(
      "Outcome",
      "Ideology",
      "Treatment",
      "Effect",
      "Std. Error",
      "P-value"
    ),
    align = c('l', 'l', 'l', 'r', 'r', 'r')
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    font_size = 11
  ) %>%
  row_spec(which(all_mfx$p.value < 0.05), bold = TRUE, color = "black") %>%
  add_header_above(c(" " = 3, "Marginal Effects" = 3)) %>%
  add_footnote(
    "Note: * p < 0.05, ** p < 0.01, *** p < 0.001",
    notation = "none"
  )

####Save Table A7
html_table_ideo
save_kable(
  html_table_ideo,
  file = here::here(thisPath(), "Figures", "Appendix_7.html"), 
  width = 8,
  height = 8
)
#Note: This code creates the html to produce appendix table 7. The table was then edited in Word to improve its appearance for the final paper.


###############################
## Interaction models for party ID (Appendix Table 8)
###############################
table(h.df$party_identity)

#Label party identity
labels_12 <- c(
  "Partido Socialista Obrero Español (PSOE)", # 1
  "Unidas Podemos", # 2
  "Podemos", # 3
  "Izquierda Unida", # 4
  "Más País", # 5
  "EQUO", # 6
  "Partido Popular (PP)", # 7
  "Vox", # 8
  "Ciudadanos", # 9
  "PACMA", #10
  "Otro", #11
  "Ninguno / NS / NC" #12  (none/don’t know/no answer)
)

#Set party levels
h.df$party_name <- factor(
  h.df$party_identity,
  levels = 1:12,
  labels = labels_12
)

h.df$pspod <- factor(
  ifelse(h.df$party_identity %in% c(1, 3), "PSOE/Podemos", "Other"),
  levels = c("Other", "PSOE/Podemos")
)

table(h.df$pspod, h.df$party_identity)

##Fit interaction models for party ID
summary(
  dem16_party.lm <- lm(
    dem16 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem2_party.lm <- lm(
    dem2 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem3_party.lm <- lm(
    dem3 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem4_party.lm <- lm(
    dem4 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem5_party.lm <- lm(
    dem5 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)
summary(
  dem7_party.lm <- lm(
    dem7 ~ treat * pspod + region + gender + age + income,
    data = h.df
  )
)

#Marginal Effects for appendix table 8
mfx_dem16_2 <- avg_slopes(dem16_party.lm, variables = "treat", by = "pspod")
mfx_dem2_2 <- avg_slopes(dem2_party.lm, variables = "treat", by = "pspod")
mfx_dem3_2 <- avg_slopes(dem3_party.lm, variables = "treat", by = "pspod")
mfx_dem4_2 <- avg_slopes(dem4_party.lm, variables = "treat", by = "pspod")
mfx_dem5_2 <- avg_slopes(dem5_party.lm, variables = "treat", by = "pspod")
mfx_dem7_2 <- avg_slopes(dem7_party.lm, variables = "treat", by = "pspod")


#combine all marginal effects into single table
all_mfx_2 <- rbind(
  mfx_dem16_2 %>% mutate(outcome = "Conceptual Democracy"),
  mfx_dem2_2 %>% mutate(outcome = "Freedom of Assembly"),
  mfx_dem3_2 %>% mutate(outcome = "Freedom of Press"),
  mfx_dem4_2 %>% mutate(outcome = "Judicial Independence"),
  mfx_dem5_2 %>% mutate(outcome = "Rule of Law"),
  mfx_dem7_2 %>% mutate(outcome = "Party Freedom")
)

#Clean up prime names for table
all_mfx_2 <- all_mfx_2 %>%
  mutate(
    contrast = case_when(
      grepl("CE Treat", contrast) ~ "Elite Corruption",
      grepl("Corr Treat", contrast) ~ "General Corruption",
      grepl("CS Treat", contrast) ~ "Systemic Corruption",
      grepl("Pol Treat", contrast) ~ "Unemployment",
      TRUE ~ contrast
    )
  )

###Clean up table
formatted_table <- all_mfx_2 %>%
  mutate(
    stars = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ ""
    ),
    estimate_display = sprintf("%.3f%s", estimate, stars),
    se_display = sprintf("(%.3f)", std.error),
    p_display = sprintf("%.3f", p.value)
  ) %>%
  select(outcome, pspod, contrast, estimate_display, se_display, p_display)

#######################################
#####TO GENERATE Appendix Table 8
#######################################
html_table_party <- formatted_table %>%
  kable(
    format = "html",
    col.names = c(
      "Outcome",
      "Incumbent",
      "Treatment",
      "Effect",
      "Std. Error",
      "P-value"
    ),
    align = c('l', 'l', 'l', 'r', 'r', 'r')
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    font_size = 11
  ) %>%
  row_spec(which(all_mfx_2$p.value < 0.05), bold = TRUE, color = "black") %>%
  add_header_above(c(" " = 3, "Marginal Effects" = 3)) %>%
  footnote(
    general = "* p < 0.05, ** p < 0.01, *** p < 0.001",
    general_title = "Note:",
    footnote_as_chunk = TRUE
  )

####Save Table A8
html_table_party
save_kable(
  html_table_party, 
  file = here::here(thisPath(), "Figures", "Appendix_8.html"),
  width = 8,
  height = 8
)

#Note: This code creates the html to produce appendix table 8. The table was then edited in Word to improve its appearance for the final paper.

